subroutine migrate(neutron, k, ke, nnt)                                            
    use define                                                             
    use util_MC         
    use ifport
    use rng   
    use wwnd                                                    
    implicit none                                                              
    include "arrays.FI"                                                        
    include "pspec.FI"                                                         
    include "mcbase.FI"  
    include "xsec.FI"                                                      
    include "CMFD_LARSEN.FI"
    real(8), parameter :: mu_min = 0.01_8                                                                               
    type(particle), intent(in) :: neutron                                      
    real(8), intent(in) :: k, ke                                             
    integer, intent(out) :: nnt                                                
    real(8) :: x, mu, r, w, wm, dist, dist_surf, alpha, trkl, test                          
    integer :: i, dir, l, lc, g, gtilde, surf, btype, comptype, comptypen, nn, type_col, isurf, ttt, inc1, inc2 
    integer :: nnw, ifmc
    real(8) :: phifmc, ftent, x0, x1, x2, mur, femt, testsum
    logical escape          
    real(8), parameter :: p=0.5, s= 0.2
    real(8) :: sigts, mud, nwt
    logical :: fly
    integer :: ret_wwnd
    integer :: icadj
    real(8), parameter :: eps = 0.01_8
    real(8) :: fx
    
    integer :: ic
    x = neutron%x                                                              
    w = neutron%w
    wm = neutron%wm                                                              
    l = neutron%l                                                              
    g = neutron%g  
    mu = neutron%mu   
#ifdef _TRPCHK
    fly = .false.
    if (l .eq. 38 ) fly = .true.
#endif
    dir = -1   
    if (mu>0.) dir = 1                                                      
    nnt = 0                                                                    
    escape = false 
    srcn(l, g) = srcn(l, g) + w
    inc1 = 0
!    call sample_direction(mu, dir)
    do while (TRUE)                                                            
        comptype = complist(l)                                                 
!        r = -1./Xtot(complist(l), g) * log(drand(0)) 
#ifdef _EXPTR
        mud = 1.
        if (l>totmesh/2) mud = -1.
        if (mud .eq. dir) then
        sigts = Xtot(complist(l), g)*(1-p*mu*mud)
        r = -1./sigts*log(getRN())
!        w = w*sigts/Xtot(complist(l), g)
        w = w *exp(-p*mu*mud*r*Xtot(complist(l), g))/(1-p*mu*mud)
        else
            r = -1./Xtot(complist(l), g) * log(getRN())                                 
        endif
!        r = -1./s/Xtot(complist(l), g) * log(getRN()) 
!        w = w * exp(-Xtot(complist(l), g)*(1-s)*r)/s                                
#else                                
        r = -1./Xtot(complist(l), g) * log(getRN())                                 
#endif
        call dist2surf(x, l, dir, surf, dist_surf)                             
        dist = r * mu * dir                                                    
        do while(dist>dist_surf)                                               
            inc2 = inc2 +1 
#ifndef _NOFMC            
            ! FMC 
            ! surface terms will be condensed in FMC routines. 
            ifmc = l+(dir+1)/2
            if (abs(mu)<mu_min) then
                phifmc =  1._8/(mu_min*0.5)*w*dir
            else
                phifmc =  1._8/mu*w*dir
            endif
            phifmc =  1._8/mu*w*dir
            phi_s(ifmc, g) = phi_s(ifmc, g)+phifmc
            mmt1st(ifmc,g) = mmt1st(ifmc,g)+phifmc*abs(mu)
            mmt2nd(ifmc,g) = mmt2nd(ifmc,g)+phifmc*(mu**2)

            ! node average terms are tallied directly. 
            lc = fi2co(l)
            x0 = x-hac_cmfd(lc-1)
            x1 = x0 + dist_surf * dir  
            mur = 1._8/mu*w*dir
            if (x0< hfmcmid(lc)) then
                if (x1<hfmcmid(lc)) then
                    phi_sa(lc,g) = phi_sa(lc,g)+g_av(lc)*mur*dist_surf
                else
                    phi_sa(lc,g) = phi_sa(lc,g)+g_av(lc)*mur*(hfmcmid(lc)-x0)
                    phi_sa(lc+1,g) = phi_sa(lc+1,g)+g_av(lc+1)*mur*(x1-hfmcmid(lc))
                endif
            else
                if (x1<hfmcmid(lc)) then
                    phi_sa(lc,g) = phi_sa(lc,g)+g_av(lc)*mur*(hfmcmid(lc)-x1)
                    phi_sa(lc+1,g) = phi_sa(lc+1,g)+g_av(lc+1)*mur*(x0-hfmcmid(lc))
                else
                    phi_sa(lc+1,g) = phi_sa(lc+1,g)+g_av(lc+1)*mur*dist_surf
                endif
            endif
            
            femt = abs( mur*( -1._8/2/h_cmfd(lc)*(x1**2-x0**2)+(x1-x0) ) )
            mmtAbsr(lc, g) = mmtAbsr(lc, g)+Xabs(comptype, g)*femt
            mmtNuF(lc, g) = mmtNuF(lc, g)+Xnuf(comptype, g)*femt
            femt = abs( mur*( 1._8/2/h_cmfd(lc)*(x1**2-x0**2) ) )
            mmtAbsr(lc+1, g) = mmtAbsr(lc+1, g)+Xabs(comptype, g)*femt
            mmtNuF(lc+1, g) = mmtNuF(lc+1, g)+Xnuf(comptype, g)*femt

    
#endif
            
            x1=x
            x = x + dist_surf * dir                                            
            x2=x
            
            if (abs(mu) < eps) then
                trkl = dist_surf*w*2/eps
            else
                trkl = dist_surf/mu*dir*w
            endif
            
            trk_est(l, g) = trk_est(l, g) + trkl  
            gktrk=gktrk+xnuf(complist(l),g)*trkl
            ! tally for gcmfd functionals
            ic=fi2co(l)

            if (cmfdtype .eq. 8) then
                call tally_gcmfd4(l, ic)
            elseif (cmfdtype .eq. 7) then
                call tally_gcmfd3(l, ic)
            elseif (cmfdtype .eq. 9) then
                call tally_gcmfd2(l,ic)
            elseif (cmfdtype .eq. 10) then
                call tally_gcmfd5(l,ic)
            endif
            
#ifdef _2NDHRMN
            trk_est2(l,g) = trk_est2(l,g) + trkl*wm/w
#endif
            phi_f2(l,g) = phi_f2(l,g) + trkl*mu**2                             
!            trk_estT(l, g) = trk_estT(l, g) + trkl                             
            isurf = (l-1)*2+(dir+1)/2+1                                       
            jf(isurf, g) = jf(isurf, g)+ w
            icadj=(dir-1)/2  ! 0: mu>0, -1:mu<0
#ifdef CMFDLARSEN
            phi1La(l,g)=phi1La(l,g)+trkl*mu
            phi2Ls(l+icadj,g)=phi2Ls(l+icadj,g)+abs(mu)*w
#endif           
    
            if (cmfdtype .eq. TGCMFD5) then ! GCMFD5
                fx=(hac_cmfd(ic)-x)/h_cmfd(ic)
                if (fx<1e-10_8) fx=0._8
                fphi2fs(1,l+icadj,g)=fphi2fs(1,l+icadj,g)+abs(mu)*w*fx
                fx=(x-hac_cmfd(ic-1))/h_cmfd(ic)
                if (fx<1e-10_8) fx=0._8
                fphi2fs(2,l+icadj,g)=fphi2fs(2,l+icadj,g)+abs(mu)*w*fx
            endif
            
            if (ic .ne. fi2co(surf)) then
                phi1c(ic+icadj,g)=phi1c(ic+icadj,g)+w*dir
            endif                                                           
            l = surf
                                                 
            dist = dist - dist_surf                                            
            if ((l.eq.0).or.(l.eq.totmesh+1)) then                             
                if (l.eq.0) then                                               
                    btype = 1 ! left                                           
                else                                                           
                    btype = 2 ! right                                          
                end if                                                         
                if (albedo(btype)<1e-10) then 
                    jf(isurf, g) = jf(isurf, g)- w    
                    if (btype.eq.1) then                                       
                        l = l+1                                                
                    else                                                       
                        l = l-1                                                
                    end if                                                     
                    mu  = mu*-1              
                    icadj=(dir-1)/2                                  
                    phi2Ls(l+icadj,g)=phi2Ls(l+icadj,g)+abs(mu)*w
                    phi1c(ic+icadj,g)=phi1c(ic+icadj,g)-w*dir
                    dir = dir * -1     
                                                            
                ! vacuum                                                       
                else                                                           
                    escape = TRUE  
                    phi1p((dir+1)/2+1,g)=phi1p((dir+1)/2+1,g)+w
                    exit                                                       
                end if                                                         
            else                                                               
#ifdef _WWND
            ret_wwnd=event_ww(w,l,g) 
            if (ret_wwnd .eq. SPLT) then
                call addntrn(ntrn(:,bat(CUR)), nq(bat(CUR)), x, l, mu, w, g)                                
            elseif (ret_wwnd .eq. TERM) then
                exit
            else
            
            endif
#endif       
                comptypen = complist(l)                                        
                if (comptypen.ne.comptype) then                                
                    comptype = comptypen                                       
!                    r = -1./Xtot(complist(l), g) * log(drand(0))
                    !r = -1./Xtot(complist(l), g) * log(GetRN())
#ifdef _EXPTR
                    mud = 1.
                    if (l>totmesh/2) mud = -1.
                    if (mud .eq. dir) then
                    sigts = Xtot(complist(l), g)*(1-p*mu*mud)
                    r = -1./sigts*log(getRN())
            !        w = w*sigts/Xtot(complist(l), g)
                    w = w *exp(-p*mu*mud*r*Xtot(complist(l), g))/(1-p*mu*mud)
                    else
                        r = -1./Xtot(complist(l), g) * log(getRN())                                 
                    endif
!                    r = -1./s/Xtot(complist(l), g) * log(getRN()) 
!                    w = w * exp(-Xtot(complist(l), g)*(1-s)*r)/s                                
#else                                
                    r = -1./Xtot(complist(l), g) * log(getRN())                                 
#endif
                    dist = r * mu * dir                                        
                end if                                                         
            end if                                                             
            call dist2surf(x, l, dir, surf, dist_surf)                         
        end do
                                                                         
        if (escape.eq.FALSE) then                                              
#ifndef _NOFMC
            ! FMC
            lc = fi2co(l)
            x0 = x-hac_cmfd(lc-1)
            x1 = x0 + dist * dir  
            mur = 1._8/mu*w*dir
            if (x0< hfmcmid(lc)) then
                if (x1<hfmcmid(lc)) then
                    phi_sa(lc,g) = phi_sa(lc,g)+g_av(lc)*mur*dist
                else
                    phi_sa(lc,g) = phi_sa(lc,g)+g_av(lc)*mur*(hfmcmid(lc)-x0)
                    phi_sa(lc+1,g) = phi_sa(lc+1,g)+g_av(lc+1)*mur*(x1-hfmcmid(lc))
                endif
            else
                if (x1<hfmcmid(lc)) then
                    phi_sa(lc,g) = phi_sa(lc,g)+g_av(lc)*mur*(hfmcmid(lc)-x1)
                    phi_sa(lc+1,g) = phi_sa(lc+1,g)+g_av(lc+1)*mur*(x0-hfmcmid(lc))
                else
                    phi_sa(lc+1,g) = phi_sa(lc+1,g)+g_av(lc+1)*mur*dist
                endif
            endif
            
            femt = abs( mur*( -1._8/2/h_cmfd(lc)*(x1**2-x0**2)+(x1-x0) ) )
            mmtAbsr(lc, g) = mmtAbsr(lc, g)+Xabs(comptype, g)*femt
            mmtNuF(lc, g) = mmtNuF(lc, g)+Xnuf(comptype, g)*femt
            femt = abs( mur*( 1._8/2/h_cmfd(lc)*(x1**2-x0**2) ) )
            mmtAbsr(lc+1, g) = mmtAbsr(lc+1, g)+Xabs(comptype, g)*femt
            mmtNuF(lc+1, g) = mmtNuF(lc+1, g)+Xnuf(comptype, g)*femt

            
#endif                                                                               
            x1=x
            x = x + dist*dir
            x2=x
!            if (dist==0. .and. mu==0.) then
!                write(*, '(a)'), ' Error'
!            end if 
            if (abs(mu) < eps) then
                trkl = dist*w*2/eps
            else
                trkl = dist/mu*dir*w
            endif
            trk_est(l, g) = trk_est(l, g) + trkl
            gktrk=gktrk+xnuf(complist(l),g)*trkl
            ic=fi2co(l)
            
!            phi0f(l,g)=phi0f(l,g)+trkl
!            phi2f(l,g)=phi2f(l,g)+trkl*mu**2 
!            fphi1_hf(1,l,g)=fphi1_hf(1,l,g)+intg_f1(ic,x1,x2)*mu*w
!            fphi1_hf(2,l,g)=fphi1_hf(2,l,g)+intg_f2(ic,x1,x2)*mu*w
!                                                               
!            gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1(ic,x1,x2)*mu**2*w  
!            gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(ic,x1,x2)*mu**2*w  
!            gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd4(ic,x1,x2)*mu**2*w  


            if (cmfdtype .eq. 8) then
                call tally_gcmfd4(l, ic)
            elseif (cmfdtype .eq. 7) then
                call tally_gcmfd3(l, ic)
            elseif (cmfdtype .eq. 9) then
                call tally_gcmfd2(l,ic)
            elseif (cmfdtype .eq. 10) then
                call tally_gcmfd5(l,ic)
            endif
            
            
#ifdef CMFDLARSEN
            phi1La(l,g)=phi1La(l,g)+trkl*mu
#endif               
                                                                                           
            
#ifdef _2NDHRMN
            trk_est2(l,g) = trk_est2(l,g) + trkl*wm/w
#endif
            phi_f2(l,g) = phi_f2(l,g) + trkl*mu**2                             
!            trk_estT(l, g) = trk_estT(l, g) + trkl
            isurf = (l-1)*2+(dir+1)/2+1                                       
            !jf(isurf, g) = jf(isurf, g)+ w                                    
            
            call sample_Nneutron(l, g, w, k, ke, nn, nnw, nwt)  
            !absn(l, g) = absn(l, g) + nn*w
            nnt = nnt+nn                                                       
            do i=1, nn                                                         
!                nqn = nqn +1                                                   
!                ntrn(nqn, batn)%x = x                                          
!                ntrn(nqn, batn)%w = 1.                                         
!                ntrn(nqn, batn)%l = l                                          
!                ntrn(nqn, batn)%g = nn_group() 
!                ntrn(nqn, batn)%mu = 1-2*drand(0)
!                call addntrn(ntrn(:,batn), nqn, x, l, 1-2*drand(0), 1._8, nn_group())    
#ifdef _2NDHRMN
#else
                call addntrn(ntrn(:,bat(NXT)), nq(bat(NXT)), x, l, 1-2* getRN() , nwt, nn_group())                                
#endif
                nfn(l)=nfn(l)+1
            end do     
            do i=1, nnw                                                        
#ifdef _2NDHRMN
#else
                call addntrn(ntrn(:,bat(CUR)), nq(bat(CUR)), x, l, 1-2* getRN() , nwt, nn_group())                                
#endif
!                call addntrn(ntrn(:,bat(CUR)), nq(bat(CUR)), x, l, 1-2* getRN() , 1._8, nn_group())                                
            end do    
            
#ifdef _TRPCHK
            if (fly .eq. .true. .and. l .eq. 88) then
                tempt(1) = tempt(1) + 1
                tempt(2) = tempt(2) + w
                if (nn>0) tempt(3) = tempt(3) + nn
            endif
#endif

#ifdef _WWND
            w = w*Xsct(comptype, g)/Xtot(comptype, g)
#else
        if (impl_cpt) then !===============================================
            w = w*Xsct(comptype, g)/Xtot(comptype, g)
            
            gtilde = sample_scattering_group(l, g)                              
            sctn(g, gtilde, l) = sctn(g, gtilde, l)+w
            g = gtilde
            ! Russian roullete

            if (w < 0.25_8) then
                if(getRN()<w/0.5_8) then
                    w = 0.5_8
                else
                    w = 0.
                    exit
                endif
            endif
        else          
            type_col = sample_collision(l, g)                                  
            if (type_col.eq.COL_AB) then    
                absn(l, g) = absn(l, g) + w                                   
                exit                                                           
            else                                                               
!                if(ng>1) then
                gtilde = sample_scattering_group(l, g)                              
                sctn(g, gtilde, l) = sctn(g, gtilde, l)+w
                g = gtilde
!                end if
            end if                                                             
        endif
        !===================================================================
#endif
        else                                                                   
            exit                                                               
        end if                                                                 
#ifdef _LANISO        
        call sample_direction_anisotropic((mu), mu, dir)
#else
        call sample_direction(mu, dir)                                         
#endif
    end do    
contains

    subroutine tally_gcmfd2(l,lc)
        integer, intent(in) :: l, lc
        integer :: typecase
        real(8) :: xs(3)
        
        phi0(lc,g)=phi0(lc,g)+trkl
        phi2(lc,g)=phi2(lc,g)+trkl*mu**2 
        call intg_f_gcmfd3(l, lc, x1, x2)
        
    endsubroutine

    subroutine tally_gcmfd3(l, lc)
    
        integer, intent(in) :: l, lc
        integer :: typecase
        real(8) :: xs(3)

        phi0(lc,g)=phi0(lc,g)+trkl
        phi2(lc,g)=phi2(lc,g)+trkl*mu**2 
        
        call intg_f_gcmfd3(l, lc, x1, x2)
        call intg_g_gcmfd3(l, lc, x1, x2)
        
    end subroutine   

    subroutine intg_f_gcmfd3(l, lc, x1, x2)
        integer, intent(in) :: l, lc
        real(8), intent(in) :: x1,x2
        
        fphi1_hf(1,l,g)=fphi1_hf(1,l,g)+intg_f1_gcmfd3(lc,x1,x2)*dir*w
        fphi1_hf(2,l,g)=fphi1_hf(2,l,g)+intg_f2_gcmfd3(lc,x1,x2)*dir*w
    endsubroutine


    subroutine intg_g_gcmfd3(l, lc, x1, x2)
        integer, intent(in) :: l, lc
        real(8), intent(in) :: x1, x2
        REAL(8) :: test1, test2
        if (lc .eq. 1) then
            gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_gLEFT_gcmfd3(lc,x1,x2)*w/abs(mu) 
            gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd3(lc,x1,x2)*w/abs(mu) 
        elseif (lc .eq. totm_cmfd) then
            gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1_gcmfd3(lc,x1,x2)*w/abs(mu) 
            gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_gRIGHT_gcmfd3(lc,x1,x2)*w/abs(mu) 
        else
            gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1_gcmfd3(lc,x1,x2)*w/abs(mu) 
            gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd3(lc,x1,x2)*w/abs(mu) 
            gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd3(lc,x1,x2)*w/abs(mu) 
        endif
    endsubroutine
    
    function intg_f1_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        ret1=1._8/2/h_cmfd(lc)*(x2-hac_cmfd(lc))**2
        ret2=1._8/2/h_cmfd(lc)*(x1-hac_cmfd(lc))**2
        ret=abs(ret1-ret2)
    end function
    function intg_f2_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        ret1=-1._8/2/h_cmfd(lc)*(x2-hac_cmfd(lc-1))**2
        ret2=-1._8/2/h_cmfd(lc)*(x1-hac_cmfd(lc-1))**2
        ret=abs(ret1-ret2)
    end function
    function intg_g1_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        ret1=1._8/(6*h_hf(lc-1)*h_cmfd(lc))*(x2-hac_cmfd(lc))**3
        ret2=1._8/(6*h_hf(lc-1)*h_cmfd(lc))*(x1-hac_cmfd(lc))**3
        ret=abs(ret1-ret2)
    end function
    function intg_g2_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        ret1=x2-1._8/6*(1._8/(h_hf(lc-1)*h_cmfd(lc))*(x2-hac_cmfd(lc))**3+ \
                       1._8/(h_hf(lc)*h_cmfd(lc))*(x2-hac_cmfd(lc-1))**3)
        ret2=x1-1._8/6*(1._8/(h_hf(lc-1)*h_cmfd(lc))*(x1-hac_cmfd(lc))**3+ \
                       1._8/(h_hf(lc)*h_cmfd(lc))*(x1-hac_cmfd(lc-1))**3)
        ret=abs(ret1-ret2)
    end function
    function intg_g3_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        ret1=1._8/(6*h_hf(lc)*h_cmfd(lc))*(x2-hac_cmfd(lc-1))**3
        ret2=1._8/(6*h_hf(lc)*h_cmfd(lc))*(x1-hac_cmfd(lc-1))**3
        ret=abs(ret1-ret2)
    end function
    
    function intg_gLEFT_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        ret1=x2-1._8/(6*h_hf(1)*h_cmfd(1))*x2**3
        ret2=x1-1._8/(6*h_hf(1)*h_cmfd(1))*x1**3
        ret=abs(ret1-ret2)
    end function

    function intg_gRIGHT_gcmfd3(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8) :: XR
        XR = hac_cmfd(totm_cmfd)
        ret1=x2+1._8/(6*h_hf(totm_cmfd-1)*h_cmfd(totm_cmfd))*(XR-x2)**3
        ret2=x1+1._8/(6*h_hf(totm_cmfd-1)*h_cmfd(totm_cmfd))*(XR-x1)**3
        ret=abs(ret1-ret2)
    end function

    
    subroutine tally_gcmfd4(l, lc)
    
        integer, intent(in) :: l, lc
        integer :: typecase
        real(8) :: xs(3)

        phi0(lc,g)=phi0(lc,g)+trkl
        phi2(lc,g)=phi2(lc,g)+trkl*mu**2 
        typecase=sort_x(lc, x1, x2, xs)
        
        call intg_f_gcmfd4(l, lc, typecase, xs)
        call intg_g_gcmfd4(l, lc, typecase, xs)
        
        !! phi2 for half mesh
        select case(typecase)
        case(1)
            phi2_hf(1,lc,g)=phi2_hf(1,lc,g)+(x2-x1)*mu*w
        case(2)
            phi2_hf(2,lc,g)=phi2_hf(2,lc,g)+(xs(3)-xs(2))*abs(mu)*w
            phi2_hf(1,lc,g)=phi2_hf(1,lc,g)+(xs(2)-xs(1))*abs(mu)*w
        case(3)
            phi2_hf(2,lc,g)=phi2_hf(2,lc,g)+(x2-x1)*mu*w
        endselect
            
    end subroutine    
    
    
    function sort_x(lc, x1, x2, xs) result(casetype)
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8), intent(out) :: xs(3)
        integer :: casetype
        real(8) :: xm
        
        xm = hcmfdhf(lc)
        if (x1<x2) then
            xs(1)=x1
            xs(2)=x2
        else
            xs(1)=x2
            xs(2)=x1
        endif
        
        if (xm>xs(1)) then
            if (xm>xs(2)) then
                casetype=1
                xs(3)=xm
            else
                casetype=2
                xs(3)=xs(2)
                xs(2)=xm
            endif
        else
            casetype=3
            xs(3)=xs(2)
            xs(2)=xs(1)
            xs(1)=xm
        endif
        
    endfunction
    
    subroutine intg_f_gcmfd4(l, lc, typecase, xs)
        integer, intent(in) :: l, lc, typecase
        real(8), intent(in) :: xs(3)
        
        select case(typecase)
        case(1)
            fphi1_hf(1,l,g)=fphi1_hf(1,l,g)+intg_f1_gcmfd4(lc,xs(1),xs(2))*dir*w
        case(2)
            fphi1_hf(1,l,g)=fphi1_hf(1,l,g)+intg_f1_gcmfd4(lc,xs(1),xs(2))*dir*w
            fphi1_hf(2,l,g)=fphi1_hf(2,l,g)+intg_f2_gcmfd4(lc,xs(2),xs(3))*dir*w
        case(3)
            fphi1_hf(2,l,g)=fphi1_hf(2,l,g)+intg_f2_gcmfd4(lc,xs(2),xs(3))*dir*w
        endselect 
    endsubroutine


    subroutine intg_g_gcmfd4(l, lc, typecase, xs)
        integer, intent(in) :: l, lc, typecase
        real(8), intent(in) :: xs(3)
        
        if (lc .eq. 1) then ! left boundary
            select case(typecase)
            case(1)
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+(xs(2)-xs(1))*w/abs(mu) 
            case(2)
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+(xs(2)-xs(1))*w/abs(mu) 
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
                gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
            case(3)
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
                gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
            endselect
        elseif (lc .eq. totm_cmfd) then ! right boundary
            select case(typecase)
            case(1)
                gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
            case(2)
                gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+(xs(3)-xs(2))*w/abs(mu) 
            case(3)
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+(xs(3)-xs(2))*w/abs(mu) 
            endselect
        else
            select case(typecase)
            case(1)
                gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
            case(2)
                gphi0_hf(1,l,g)=gphi0_hf(1,l,g)+intg_g1_gcmfd4(lc,xs(1),xs(2))*w/abs(mu) 
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(1),xs(3))*w/abs(mu) 
                gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
            case(3)
                gphi0_hf(2,l,g)=gphi0_hf(2,l,g)+intg_g2_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
                gphi0_hf(3,l,g)=gphi0_hf(3,l,g)+intg_g3_gcmfd4(lc,xs(2),xs(3))*w/abs(mu) 
            endselect
        endif
    
    endsubroutine

    function intg_f1_gcmfd4(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8) :: xm
        xm = hcmfdhf(lc)
        ret1=-1._8/h_cmfd(lc)*(x2-xm)**2
        ret2=-1._8/h_cmfd(lc)*(x1-xm)**2
        ret=ret1-ret2
    end function
    function intg_f2_gcmfd4(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8) :: xm
        xm = hcmfdhf(lc)
        ret1=1._8/h_cmfd(lc)*(x2-xm)**2
        ret2=1._8/h_cmfd(lc)*(x1-xm)**2
        ret=ret1-ret2
    end function
    function intg_g1_gcmfd4(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8) :: xm
        xm = hcmfdhf(lc)
        ret1=2._8/(3*h_cmfd(lc)**2)*(x1-xm)**3
        ret2=2._8/(3*h_cmfd(lc)**2)*(x2-xm)**3
        ret=abs(ret1-ret2)
    end function
    function intg_g2_gcmfd4(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8) :: xm
        xm = hcmfdhf(lc)
        ret1=x1-2._8/(3*h_cmfd(lc)**2)*(x1-xm)**3
        ret2=x2-2._8/(3*h_cmfd(lc)**2)*(x2-xm)**3
        ret=abs(ret1-ret2)
    end function
    function intg_g3_gcmfd4(lc,x1,x2) result(ret)
        real(8) :: ret, ret1, ret2
        integer, intent(in) :: lc
        real(8), intent(in) :: x1, x2
        real(8) :: xm
        xm = hcmfdhf(lc)
        ret1=2._8/(3*h_cmfd(lc)**2)*(x1-xm)**3
        ret2=2._8/(3*h_cmfd(lc)**2)*(x2-xm)**3
        ret=abs(ret1-ret2)
    end function
    
    
    subroutine tally_gcmfd5(l, lc)
    
        integer, intent(in) :: l, lc
        integer :: typecase
        real(8) :: xs(3)

        phi2f(l,g)=phi2f(l,g)+trkl*mu**2 
        phi2(lc,g)=phi2(lc,g)+trkl*mu**2 
        phi0(lc,g)=phi0(lc,g)+trkl
        
        call intg_f_gcmfd3(l, lc, x1, x2)
        call intg_g_gcmfd3(l, lc, x1, x2)
        
    end subroutine   

    
end subroutine                                                                 
                                                                               